COMMENT ⊗   VALID 00004 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	begin comment tests on brownian motion
C00004 00003	here is a recursive procedure to find where a random brownian motion
C00008 00004	this is the driver program
C00010 ENDMK
C⊗;
begin comment tests on brownian motion;
external procedure bail;
require 2000 system_pdl;

comment the first subroutines are for generating normal deviates;

real nrm; comment if nonzero, this is the next normal deviate;

real procedure normdev;
begin real x,v1,v2,s; comment this is algorithm 3.4.1P;
if nrm then
	begin x←nrm; nrm←0.0; return(x);
	end;
while true do
	begin v1←2*ran(0)-1; v2←2*ran(0)-1; s←v1↑2+v2↑2;
	if s<1 then
		begin x←sqrt(-2*log(s)/s);
		nrm←v1*x; return(v2*x);
		end;
	end;
end;
comment here is a recursive procedure to find where a random brownian motion
crosses the line b(x)=0;

integer array count,zeros[0:25]; comment instrumentation to measure efficiency;
integer unicount,normcount; comment the number of random numbers generated;

real array delta,scale[0:25]; comment delta[k]=1/2↑{k+1}, scale[k]=1/2↑{k/2+1};

recursive boolean procedure Bzero(real y,a,b,r,s; integer k,m);
begin comment this procedure locates all zeroes between y and y+2↑-k, with
resolution up to 2↑m, for a Brownian function A(x) such that A(y)=a and
A(y+2↑-k)=b. There is also a global function B(x) such that we have
either B(x)=sA(x) or B(x)=reflection(sA(x)) in the interval y≤x≤y+2↑-k,
with the latter (reflection) alternative holding if and only if r=-1 and there
is at least one zero in the interval. (We have |r|=|s|=1.) The procedure returns
the value false if it finds no zeroes, otherwise it returns true;
real c;
count[k]←count[k]+1;
if a*b>0 then
	begin unicount←unicount+1;
	if ran(0)<exp(-a*b/delta[k]) then
		begin b←-b; r←-r;
		end
	else return(false);
	end;
comment print('15&'12,"level ",k,": B[",cvf(y),"]=",cvf(s*a),
		", B[",cvf(y+2*delta[k]),"]=",cvf(r*s*b));
if r*a*b>0 then zeros[k]←zeros[k]+2 else zeros[k]←zeros[k]+1;
if k<m then
	begin c←(a+b)/2+scale[k]*normdev; normcount←normcount+1;
	if Bzero(y,a,c,r,s,k+1,m) then Bzero(y+delta[k],c,b,1,r*s,k+1,m)
	else Bzero(y+delta[k],c,b,r,s,k+1,m);
	end;
return(true);
end;

procedure Brownianzeros(real a; integer m,s);
begin comment this procedure initializes memory and calls the recursive
Bzero procedure to find all zeroes of a Brownian function B(x) with B(0)=a,
to resolution 2↑m. The random number generator starts with seed value s;
integer k; 
unicount←0; normcount←1;
for k←0 step 1 until m do 
	begin count[k]←zeros[k]←0;
	delta[k]←2↑(-1-k); scale[k]←sqrt(delta[k]/2);
	end;
ran(s); nrm←0.0;
Bzero(0,a,a+normdev,1,1,0,m);
end;
comment this is the driver program;

integer array gcount,gzeros[0:25]; comment grand totals;
integer gunicount,gnormcount; comment likewise;
integer j,k,m,reps,seedoffset;

gunicount←gnormcount←0; arrclr(gcount); arrclr(gzeros);
setprint("brown.tmp","b");
m←20; reps←100; seedoffset←1000; bail;

for j←1 step 1 until reps do
	begin print('15&'12,"Beginning of experiment ",j);
	Brownianzeros(2/3,m,j+seedoffset);
	print('15&'12,"Number of iterations and number of zeros:");
	for k←1 step 1 until m do
		begin print('15&'12,k,":   ",count[k],"   ",zeros[k]);
		gcount[k]←gcount[k]+count[k];
		gzeros[k]←gzeros[k]+zeros[k];
		end;
	print('15&'12,"Random deviates: ",unicount," uniform, ",
		normcount," normal.");
	gunicount←gunicount+unicount;
	gnormcount←gnormcount+normcount;
	end;
print('15&'12,"Grand totals:");
for k←1 step 1 until m do
	print('15&'12,k,":   ",gcount[k],"   ",gzeros[k]);
print('15&'12,"Random deviates: ",gunicount," uniform, ",
	gnormcount," normal.");
end